%===============Boyi Li 20221104================%
%=========MY==ValoMC+Model+K-Wave+2DTR==========%
%======For calculaton absorbed+K-Wave+2DTR======%

% Compute the absorbed optical energy density.
% multiply by
% 1e6   to convert [J/mm^2] to [J/m^2]
% 1e-3  to set the total energy in the pulse to 1 mJ
%


PML_size=10;

vmcmedium.absorbed_energy = vmcmedium.absorption_coefficient .* solution.grid_fluence*1e3; % [J/m3]

% Compute the initial pressure distribution
source.p0 = vmcmedium.gruneisen_parameter .* vmcmedium.absorbed_energy;  % [Pa]



sensor.mask = zeros(kgrid.Nx, kgrid.Ny);
sensor.mask(118:182,10) = 1;


[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed/2);
input_args = {'PMLInside', false, 'PMLSize', PML_size, 'Smooth', false, 'PlotPML', false};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, 'PMLInside', false);
ylabel('Axial-Dirction [mm]');
xlabel('Lateral-Dirction [mm]');

figure('rend','painters','pos',[10 10 1200 400])

% plot the initial pressure and sensor distribution
subplot(1,3,1)
hold on

% We have to swap x and y again
imagesc( kgrid.x_vec*1e3, kgrid.y_vec*1e3, rot90(vmcmedium.absorbed_energy), [min(vmcmedium.absorbed_energy(:)) ...
        max(vmcmedium.absorbed_energy(:))]);

ylabel('Axial-Dirction [mm]');
xlabel('Lateral-Dirction [mm]');

colormap default;

c = colorbar;  % create a colorbar
axis image;
title('Absorbed energy [J/m3]');
hold off

subplot(1,3,2)

imagesc(kgrid.x_vec*1e3, kgrid.y_vec*1e3, transpose(source.p0 + ...
        cart2grid(kgrid, sensor.mask)), [min(source.p0(:)) ...
        max(source.p0(:))]);


colormap(getColorMap);
ylabel('Axial-Dirction [mm]');
xlabel('Lateral-Dirction [mm]');
c = colorbar;  % create a colorbar
axis image;
title('Initial pressure [Pa]');
hold off


sensor.time_reversal_boundary_data = sensor_data;

% run the time reversal reconstruction
p0_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor);
subplot(1,3,3)
imagesc(p0_recon')
colormap(getColorMap);

c = colorbar;  % create a colorbar
axis image;
title('REC');
ylabel('Axial-Dirction [mm]');
xlabel('Lateral-Dirction [mm]');
hold off



% plot the simulated sensor data
figure,
imagesc(flip(rot90(sensor_data)),[-0.05,0.25]);
colormap(getColorMap);
xlabel('Sensor Position');
ylabel('Time Step');
c = colorbar;  % create a colorbar
colorbar;
title('Sensor data');



